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Abstract. In this paper, we study the accuracy of Irving-Kirkwood type of formulas 
for the approximation of continuum quantities from atomistic simulations. Such 
formulas are derived by expressing the displacement, deformation gradient and stress 
in terms of certain kernel functions. We propose two criteria for choosing the 
kernel functions to significantly improve the sampling accuracy. We present a simple 
procedure to construct kernel functions that meet these criteria. Further, numerical 
tests on homogeneous and non-homogeneous systems provide validations for our 
analysis. 
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1. Introduction 

Molecular models have undoubtedly been one of the most important methods for 
material simulations. By following particle trajectories, molecular simulations generate 
a large amount of data that need to be further processed to extract averaged 
quantities. For solid systems, of particular interest are quantities that correspond 
to variables in elasticity models, e.g., displacement, strain and stress. Therefore, 
it is of great practical interest to develop appropriate formulas that map particle 
trajectories to continuum variables while still following fundamental principles in 
continuum mechanics. For homogeneous systems, the stress can be obtained from the 
second law of thermodynamics [1, 2, 3, 4, 5, 6, 7, 8]. The stress obtained this way is 
known as the virial stress. For non-homogeneous systems, the local stress (in space and 
time) can be defined following the seminal work of Irving and Kirkwood [9] , where point 
distributions were introduced to represent a continuous displacement and velocity field. 
The formula for the stress is derived based on conservation laws, which is consistent 
with the continuum mechanics framework. A more practical approach was proposed by 
Hardy [10], where a smooth kernel function is used to replace the point distributions. 
This approach leads to a spatially averaged stress, called the Hardy stress, derived 
also from conservation laws. Similar work can be found in [11, 12, 13, 14, 15, 16, 17]. 
Recently, Yang et al [18] proposed to further generalized the formulas to involve the 
temporal and spatial averages simultaneously. 

Zimmerman [19] discussed several practical issues with the calculation of the Hardy 
stress, including the choice of kernel functions and the sample radius. However, to our 
knowledge, there has been no quantitative estimate of the sampling error. In this paper, 
we analyze the spatial sampling error at zero temperature and relate the error to certain 
moment conditions for the kernel function. Further, we suggest a simple procedure to 
construct hybrid kernel functions to satisfy the moment conditions and to significantly 
improve the sampling accuracy. 

To study the accuracy within a standard numerical analysis framework, we have in 
mind a smooth function, with values only given at atomic positions. Our goal is then to 
reconstruct the smooth function with maximal order of accuracy. We employ a kernel 
function to form a set of basis functions centered at each atom, and use the values of 
the function as the nodal values. Using Taylor expansions, we have obtained the leading 
terms in the error. By choosing kernel functions for which the leading terms are small 
or zero, we maximize the order of accuracy. We will show that these kernel functions 
also lead to more accurate formulas for the stress calculation. 

The rest of the paper is organized as follows. In section 2, we show a simple 
derivation of the Hardy stress for molecular mechanics models. In section 3, we provide 
error estimate for the displacement, deformation gradient and the stress. Then in section 
4, a procedure for constructing accurate kernel functions is presented. We show some 
numerical tests in section 5. 
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2. The Irving-Kirkwood Type of Formulas 

In this section, we briefly describe how continuum quantities can be defined from the 
atomic positions. We will focus primarily on the formulas proposed by Hardy [10]. These 
formulas can be derived from a molecular dynamics model based on the conservation 
of momentum and energy [10]. Here, we present a similar derivation starting with a 
molecular mechanics model, given by 

-W x V + b t = 0, (1) 

where bi is a body force. We will let fi = — S7 Xi V be the inter-molecular force on atom 
i, and we assume that fi admits the following decomposition, 

Ji = / j Tij-i Jij = ~Jji- {^) 

This leads to 

- E ^ = **■ ( 3 ) 

Now we define a continuously distributed force, 

b{x) = ^b jV {x- Xj ). (4) 

3 

This formula provides a means to sample the body force at any point in space. The 
kernel function is typically assumed to be symmetric with compact support, i.e., 

ip(x) = <p(-x), (5) 

and 

/ ip(x)dx = 1. (6) 

In practice, we may start with a non-dimensional function tp (x), and (po{ x ) — when 
\x\ > 1. Then we choose a sample radius R c , and let 

1 x 

Multiplying the equation (3) by (p(x — Xj), we find that, 

-V-a = b(x), (7) 

where the expression cr is given by, 

*W = ~2 S Yl hi ® x ij b ij( x )i ( 8 ) 

and, 

bij(x) = ¥ (x - (Xi - \Xij)) d\. (9) 

Jo 

Here Xij = x,i — Xj. The tensor cr is identified as the first Piola-Kirchhoff stress since 

equation (7) has the same form as the continuum elastostatics model. The atomic 
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expression of the stress (8) is in the same form as the Hardy stress derived from molecular 
dynamics models [10, 21, 22]. 

However, this formulation only involves body forces and stress. The formula for 
the displacement is not part of this formulation. One may use a formula similar to (4) 
to define the local displacement. However, it turns out that the formula is inconsistent 
with the continuum elastodynamics model. In this case, we have to rely on the Irving- 
Kirkwood formalism, where the local velocity is given by, 

v(x, t) = - ^2 rriiVi(t)(p(x - x { ), (10) 

" i 

in which, 

p(x) = ^mjipjx - Xj). (11) 

i 

In order to be consistent with the definition of local velocity, we define the 
displacement as follows, 

u(x, t) = - ^2 miUi(t)(p(x - Xi). (12) 

i 

d „ 
As a result, we have, — u = v. In addition, by taking the time derivative of (10), we 

at 
obtain an equation that resembles the elastodynamics equation, 

(13) 















! a is also given by (7). 
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q(x) = y~] miUiip 


(x-xi). 











(14) 



_ (15) 

i 

In this work, we are interested in the case when the system is made of a single 
crystal, for which TOj = m, and the underlying lattice is a simple lattice. This is the 
case where the accuracy can be significantly improved. 

3. Error Estimate 

We now turn to the error estimate of the formulas presented in the previous section. 
These estimates will be summarized as several theorems. These theorems are presented 
not for the purpose of mathematical rigor, but for the purpose of clarity. 
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3.1. Approximation of a function at the lattice points 

We notice that both of the formulas (11) and (15) are in the following form, 

p w(x) = NJ Wj<p(x — Xj). (16) 

3 

Here the function tp can be considered as a regularization of the delta function to 
smoothly sample data over a set of discrete points, po is the density in the reference 
lattice, po — ^o *> with ^o being the volume of the unit cell. We will always assume 
that (p is C 1 with compact support, and it satisfies (5) and (6). 

To formulate the problem as an interpolation problem in standard numerical 
analysis, we suppose that lUj corresponds to a smooth function, for instance, Wi = w(xi), 
where w(x), typically unknown, has enough smoothness. In addition, the lattice 
spacing, denoted by e, is regarded as a small parameter. It means that over the atomic 
scale, w is a slowly varying function. In this case, we wish to be able to reproduce 
w(xi). Namely w(x) ~ w(x) with high accuracy. To this end, we first expand w(x) 
in Taylor series, 

w(xj) = w(Xi) + (Xj — Xi) ■ Vltf(sCj) 

1 (17) 

Combining the equations (16) and (17) and collecting terms, we found 

p w(Xi) = ^(fiXi - Xj)w(Xi) 

I (18) 

+ -^ {Xl - Xj){{Xl - X] ).Vfw{^) + 0{e% 

j 

In such an expansion, we can eliminate terms of the following form, 

^2if(Xi -Xj)(Xi- Xj ) n , 
j 
for odd powers of n. This is due to the symmetry of the lattice and the kernel function 
ip. 

Using the translational symmetry, we may define the moments, 

m o = Yl ^fo')' m % P = Yl x< j xl j^( x 3')- (I 9 ) 

3 j 

Here x^ denotes a component of the coordinate of the jth atom. As a result we have, 

p w(xi) = m w(xi) + -m 2 ■ V 2 w(x,j) + 0(e 4 ). (20) 

This leads to the following error estimate. 
Theorem 3.1. Suppose that w G C 2 , and 

m = po, (21) 

then, 

w{ Xi ) = w{ Xl ) + 0{t 2 ). (22) 
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If we further assume that w E C 4 and that 

m 2 = 0, (23) 

then, 

w(x i ) = w(x l ) + 0(e 4 ). (24) 

The conditions (21) and (23) will be referred to as moment conditions. Although 
the first moment condition (21) can be viewed as a quadrature approximation of the 
equation (6), in general (6) does not imply that the first moment condition is exactly 
satisfied. This has been observed in [23, 24]. In practice, when the sample radius R c is 
small, itlq may not be close to po, and the lemma indicates that large error may occur. 

3.2. Approximation of the gradient 

Taking the gradient of (16), we find that, 

p Vw(x) = 2_ J w j'^ i P( x ~ x j)- (25) 

j 

Inserting (17) into the equation above, we find that, 

p Vw(xi) = w(Xi) ^2 ^<-P( x i ~ x j) 

j _ (26) 

+ 2J V^(ccj - Xj) ® (xj - Xi)Vw(xi) + C(e 3 ). 

3 

We can eliminate the first term using the antisymmetry of Vy?. Using the translational 
symmetry, we define another moment, 

^i = ~ 5Z dx a( P( x j) x j- (27) 

3 

This leads to: 

Lemma 3.2. Suppose that w G C 3 (1R) 7 and 

Hi = Pol, (28) 

where I is the identity matrix, then 

X7w(xi) = Vw(Xi) + C(e 3 ). (29) 

Enforcing (28) in addition to the moment conditions (21) and (23) would certainly 
impose too many constraints on the kernel function. Notice that p\ corresponds to the 

integral / V</?(a?) ®xdx. In fact, the condition (28) is a discrete analog of the identity, 

V(p(x) ® xdx = -I, (30) 

which can be derived from (6). Therefore the moment conditions (21) and (28) are not 
independent constraints. 
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To illustrate the connection between the moment conditions (21) and (28), we 
express each lattice point in terms of three basis vectors, 6i, b 2 and b 3 , 

x = xibi + x 2 b 2 + x 3 b 3 . (31) 

Here x\ = ie, x 2 = je, and x 3 = ke, with i,j, and k being integers. We also need the 
reciprocal basis &, with the orthogonality condition, 

sm ■ o n = o m ^ n , m, n = 1, Zj 6. \3Z) 

Next we let ip(xi,x 2 ,x 3 ) = tp(x). A direct calculation yields, 

Vip(x) = d xi ^! + d X2 ip£ 2 + d X3 ip£ 3 . (33) 

Therefore, the right hand side of (27) will contain terms in the form of, 

Y Y Y d ^( x u x 2 , x 3 )x n . 

i j k 

For example, we can estimate, 

Y Y Y dx^fa, x 2 , x 3 )xi 

i j k 

=e Y Y Y ^i^( ie ' ^ e ' ke ) 1 

j k i 

=YYY^(( i + 1 } e >i e > ke }~^( ie >J e > ke ^ i (34) 

j k i 



\^ 2 YYY d ¥^^keY + 0^). 

j k i 



In the last step, we have applied a summation-by-parts along the 61 direction, and then 
switched the summation back to the summation over x. Further, the second term has 
been eliminated using the symmetry of the lattice and the function tp. Following the 
same procedure, one can show that 

Y Y Y d *M x ^ x 2, x s)x 2 = C(e 3 ). 

i j k 

Combining these results, we find that, 

/xi = -m I + C(e 3 ). (35) 

As a result, we have, 

Theorem 3.3. If the condition (28) in Lemma 3.2 is replaced by the moment condition 
(21), the same results hold. 
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3. 3. Approximation of the stress 

In this section, we focus on the implication of the moment conditions to the accuracy 
of the calculated stress. In order to assess the accuracy, we need to define the exact 
stress. This is a concept subject to a great deal of debate [25]. However, for a uniformly 
deformed system, with either infinite size or periodic boundary conditions, the stress is 
well defined, and it is given by the formula, 

3 

Here Qq is the volume of a primitive cell. Since the system is uniform, the stress does 
not depend on the choice of %. We have, 

Theorem 3.4. For a system with uniform deformation gradient, one has the error 
estimate: 

&(x i ) = a + 0(e 2 ), (37) 

under the condition (21). 

Proof. We begin with the expression of the stress (8) and (9). For a system with 
uniform deformation gradient, it can be easily verified that the inversion symmetry is 
still preserved. In addition, the force fy only depends on the relative position of the 
two atoms. Therefore, it is enough to consider the forces around a particular atom, say, 
0. Then we may rewrite the stress as, 

1 f 1 

tr{x) = -- ^ foj <H> Xoj ^2 / f{x -Xi + \x 0j )d\. (38) 

j i ° 

Using the Taylor expansion, 

if(x — Xi + AjEqj) = if(x — Xi) + \x 0j ■ V<p(x — X,j) 

+ -\ 2 (x 0j ■ V(p) 2 (x -Xi)^ , 



&( x ) = —^y] foj ® x 0j ^2(p(x-xi 



3 



we get, 

rr(x) = — 
2 

~ I 5Z f°i ® x °i 5Z X °J ' ^V( x ~ x i) + C(e 2 )- 

3 i 

If x is a lattice point, and the moment condition (21) is satisfied, then, 

~2 Yl fa ® x v Yl v( x ~ Xi ^ = ~^" J2 fa x °i = a 

3 i 3 

The second term on the right hand side becomes zero because of the inversion symmetry. 
This completes the proof. □ 
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4. Constructing Hybrid Kernels with High Order of Accuracy 

Usually, the kernel functions suggested by previous works [10, 22, 18] do not satisfy 
the moment conditions exactly. Therefore, the accuracy is not guaranteed. To ensure 
accuracy, we will construct hybrid kernel functions from existing ones. Notice that the 
second moment condition (23) has 9 equations. However, for cubic crystalline systems, 
these equations can be significantly reduced, 

Theorem 4.1. For cubic systems (b.c.c and f.c.c) the second moment is CI, where C 
is a constant and I is the 3x3 identity matrix. 

This can be proved by expressing the coordinate of each atom in terms of the basis 
vectors, and use the symmetry of the cubic lattice and the kernel function. For brevity, 
we will skip the proof. 

We now describe the procedure to construct hybrid kernel functions. To begin with, 
let ipi and (p 2 be two kernel functions that satisfy (6) and the symmetry (5). We define, 

^hybrid = A ± Lpi + A 2 (fi 2 . (39) 

with the two constants A\ and A 2 satisfying, 

A l + A 2 = 1, (40) 

in order for y?h y brid to satisfy (6). Further, to satisfy the moment condition (21) exactly, 
we need, 

A 1 ^2(p 1 (x i )+A 2 ^2<f 2 (xi) =po. (41) 

i i 

These two coupled equations can be solved to obtain the two parameters. A similar 
procedure can be applied to incorporate the second moment condition (23). 

For the purpose of this construction, it is important to have a large collection of 
independent kernel functions to work with. Some examples are provided here. They are 
already scaled properly to satisfy (6). 

(i) Cubic spline function. 

&.)=[ Sl 1 + (2 ''- 3 >'' 2 l' '' £1 - (42) 

0, r > 1. 



Here r = \ x\ + x\ + x\ is the length of the vector. 

(ii) A regularized step function [18]. 

1 0.1 

<f%(x)=\ ^ exp ^3I> r ^' (43) 

0, r > 1. 

Here r = \l x\ + x\ + x§, c ~ 2.77442 is a normalization constant. 
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(iii) Cosine kernels [22] 

' 1 3 

- J^[(l + cos7rei), \xi\ < l,i = 1,2,3, 



^(x) 



i=l 



o, 



(iv) Gaussian kernel [10] 



<P?(x) 



^0.997V27r 



otherwise. 



„ 9x?+9xZ+9x% 

f(e a ), |xi| < 1,2 = 1,2,3, 



0. 



otherwise. 



A truncated polynomial function [21] 



Vcfte) 



15, 
16' 



Y[(l-2x 2 i +xt), \xi\ < l,i = l,2,3, 



i=i 



otherwise. 
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(44) 



(45) 



(46) 



5. Numerical Experiments 

We first consider a displacement field given by the following elementary functions, 
defined on a f.c.c. lattice, with lattice constant e, 

2e 



9o{x 



9i{x 



92{X 



93{X 



9a{x 



95{X 



100' 
2e ,2x. 



100 
2e 

loo 1 

2e 

loo 1 

2e 

loo 1 

2e 



e 
2au, 

~e' 
2x, 

e 
2x, 



,2nx , 



sin 



100 e 

For each of these cases, we compute the approximation function w(x) using all the five 
kernel functions. The cut-off distance of the kernel function is set to be R c = 4e and 
the simulation box is of the size 20e x 20e x lOe. The error, measured in the maximum 
norm, is listed in table 1. We observe that starting from #2, the approximation error 
becomes large. We have chosen R c large enough so that the first moment condition (21) 
is satisfied. Also shown in the table is the second moment of each kernel. One can see 
that the moment is away from zero. 

We now construct a hybrid kernel function from cp 1 and Lp n . For an iron- alpha 
system with b.c.c. structure and lattice spacing e = 2.865(A), we have listed the results 
in table 2. For an aluminum system with f.c.c structure and lattice spacing e = 4.032(A), 
the results are listed in table 3. Clearly, from the first two rows, we can see that the 
error is nearly zero when the first moment condition (21) is satisfied, confirming that 
each kernel function will exactly reproduce constant and linear functions. From the 
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Table 1. Error ||i0(a;) — w(x)||oo without the moment conditions (23). 



w 



w 



w z \ 



w 



w n \ 



w 



w IH \ 



w 



W IV \ 



\w 



w 



go 


0.000000E+00 


0.000000E+00 


0.000000E+00 


0.000000E+00 


0.000000E+00 


9i 


0.288384E-11 


0.352140E-11 


0.313907E-11 


0.288384E-11 


0.243145E-11 


92 


0.192054E+04 


0.191935E+04 


0.186910E+04 


0.192054E+04 


0.188184E+04 


9^ 


0.134305E+06 


0.134221E+06 


0.130707E+06 


0.134305E+06 


0.131599E+06 


9^ 


0.842426E+07 


0.841901E+07 


0.819845E+07 


0.842426E+07 


0.825439E+07 


9i 


0.265355E+03 


0.265257E+03 


0.261305E+03 


0.265355E+03 


0.262241E+03 


m 2 


0.605382E+01 


0.604951E+01 


0.586577E+01 


0.605382E+01 


0.591210E+01 



third and forth rows in tables 1, 2 and 3, we see that with the second moment condition 
(23), we are also able to approximate accurately quadratic and cubic functions. 

Table 2. Displacement error for a b.c.c. system. 



w 



\W - W Hybrid \\oc 



\W - u/lloo 



\w - tu /// || 00 



go 


0.10824674E-14 


0.11102230E-15 


0.27755576E-16 


9i 


0.52735594E-15 


0.16653345E-15 


0.13877788E-15 


92 


0.47184479E-15 


0.13825538E-01 


0.16859664E-01 


9?, 


0.47184479E-15 


0.13825538E-01 


0.16859664E-01 


9i 


0.23260274E-02 


0.85314910E-01 


0.10454843E+00 


95 


0.56716289E-03 


0.16469294E-01 


0.19959146E-01 



Table 3. Displacement error for a f.c.c. system. 



W \\W - W Hybrid \\oo 



\w — it> 7 | 



\w-w in \ 



go 

9i 

92 

9?, 
9i 

95 



0.15265567E-14 
0.74940054E-15 
0.63837824E-15 
0.47184479E-15 
0.23340091E-02 
0.56909837E-03 



0.99920072E-15 
0.11379786E-14 
0.13831333E-01 
0.41493998E-01 
0.85350211E-01 
0.16476307E-01 



0.10269563E-14 

0.49960036E-15 

0.16859664E-01 

0.50578993E-01 

0.10454843E+00 

0.19959146E-01 



In the next numerical experiment, we focus only on the first moment condition 
(21). We construct hybrid kernel functions mentioned above. We solve the linear system 
composed of (40) and (41). In the case when the moments of tpi and y?2 are close, a 
combination may not satisfy the condition (21). In this case, we choose the hybrid 
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8.05 



7.95 



7.9 



7.85 



7.8 



7.75 



7.7 




15 

R C (A) 



30 



Figure 1. The first moment of tp , (p and the hybrid kernel based on ip and 
ip versus the sampling radius (R c )- 



35.5 



S 34.5 



33.5 




Figure 2. Stress for ip , ip and the hybrid kernel based on ip 
the sampling radius (R c ) with 1% stretch in b.c.c Fe system 



/// 



and ip versus 



kernel to be the one whose moment is closer to po, see Figure 1. Using this hybrid 
kernel, we compute the stress of a b.c.c Fe system under uniform 1% stretch. The 
results are displayed in Figure 2. Clearly the hybrid kernel gives much more accurate 
results, especially when the sample radius R c > 6.5A. More importantly, the error in 
this regime is quite uniformly controlled, while the error from a single kernel function 
suffers from significant fluctuations. 
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To study the accuracy of the computed strain, we computed the approximate 
gradient and made comparisons with the exact values. The results for the b.c.c system 
is shown in table 4, and for the f.c.c. system, the results are summarized in table 5. In 
both cases, the hybrid kernel is constructed from the first and third kernels, we see that 
the hybrid approach has significantly reduced the approximation error. 

Table 4. Deformation gradient error for a b.c.c system. 



w 


|| Vw - Viv Hvbrid \\ 0O 


llVto-V^Hoo 


WVw-Vw^W 


9o 


0.40661405E-16 


0.28940627E-17 


0.17948354E-17 


9i 


0.44929338E-15 


0.49453648E-02 


0.50644210E-02 


92 


0.92287289E-15 


0.98907297E-02 


0.10128842E-01 


93 


0.59230612E-02 


0.13790267E-01 


0.13979664E-01 


9a 


0.23692245E-01 


0. 15598 148E-01 


0.15403288E-01 


95 


0.84147098E-02 


0.84147098E-02 


0.84147098E-02 



Table 5. Error of the deformation gradient for a f.c.c system. 



w 


|| Vw - Viv Hybrid \\ 00 


llVw-Vt^Hoo 


||Vw- Vu/^Hoo 


9o 


0.27415087E-16 


0.47228050E-17 


0.24982699E-17 


9\ 


0.36429193E-16 


0.64767543E-05 


0.12884199E-03 


92 


0.34694470E-16 


0.12953509E-04 


0.25768398E-03 


9s 


0.20740369E-02 


0.20757141E-02 


0.20406720E-02 


9\ 


0.82961474E-02 


0.82510422E-02 


0.91934239E-02 


95 


0.84147098E-02 


0.84147098E-02 


0.84147098E-02 



Finally, we consider a crack in a b.c.c. crystal of iron-alpha. The system 
studied consists of a 3D rectangular sample, with three orthogonal axis being along 
[110], [110] and [001] directions respectively. The crack is oriented along the [001] [110] 
directions. The crack is initialized using the anisotropic linear elasticity solution [26]. 
This displacement field from the analytical expressions will be regarded as the exact 
solution. The interatomic potential is the embedded atom potential [27]. Based on 
the displacement of the atoms, we have computed the averaged displacement at each 
atomic position, and the error is listed in the table 6. Notice that in general, the averaged 
displacement given by the Hardy's formulas does not agree with the actual displacement 
at that point, as discussed in section 3. The error is large when the displacement field is 
not smooth, as confirmed by Figure 5, in which large error is observed near the crack-tip. 
In contrast, the hybrid kernel function yields much more accurate results. Finally, in 
Figure 5, we show the stress computed using the original kernels and the hybrid kernel. 
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Table 6. The comparison among ip , ip and the hybrid kernel, 
measured by different norms around a crack in a b.c.c. iron system. 



Norm 



W 



H> 



Jii 



Hybrid ip 1 and cp 



in 



|i 2.306150E+00 
| 2 5.257212E-02 
„o 2.840000E-03 



2.806700E+00 
6.393728E-02 
3.380000E-03 



2.640000E-02 
2.032535E-03 
3.700000E-04 






Figure 3. The error of the displacement (112) for the crack problem. From top to 
bottom: results computed from </?i, <^2, and the hybrid kernel. 
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Figure 4. The sampled stress <7n in the crack problem. From top to bottom: ipi, if2, 
and the hybrid kernel. 
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